pkg_list = c('zoo', 'tseries', 'fGarch','PEIP','tidyverse','gridExtra', 'gdata', 'xtable','vioplot')
# ensure existing required packages are up to date:
#update.packages(ask=FALSE, oldPkgs=pkg_list)
# Install packages if needed
for (pkg in pkg_list)
{
# Try loading the library.
if ( ! library(pkg, logical.return=TRUE, character.only=TRUE) )
{
# If the library cannot be loaded, install it; then load.
install.packages(pkg, repos = "http://cran.us.r-project.org")
library(pkg, character.only=TRUE)
}
}
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Loading required package: timeDate
## Loading required package: timeSeries
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
## Loading required package: fBasics
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 1.0.0
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks timeSeries::filter(), stats::filter()
## x dplyr::lag() masks timeSeries::lag(), stats::lag()
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
##
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
##
## Attaching package: 'gdata'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:dplyr':
##
## combine, first, last
## The following object is masked from 'package:purrr':
##
## keep
## The following object is masked from 'package:stats':
##
## nobs
## The following object is masked from 'package:utils':
##
## object.size
## The following object is masked from 'package:base':
##
## startsWith
##
## Attaching package: 'xtable'
## The following object is masked from 'package:timeSeries':
##
## align
## The following object is masked from 'package:timeDate':
##
## align
## Installing package into '/Users/hoqueme/Library/R/4.0/library'
## (as 'lib' is unspecified)
## also installing the dependency 'sm'
##
## The downloaded binary packages are in
## /var/folders/j9/dfgsrh6n401466rqzvbmlmq80000gn/T//RtmpYIN9pl/downloaded_packages
## Loading required package: sm
## Package 'sm', version 2.2-5.6: type help(sm) for summary information
#library(fBasics)
rm(list=ls())# Remove objects from enviornment
dateStart <- "2005-01-01"
dateEnd <- "2018-11-28"
# Apple data
APPLE <- get.hist.quote(instrument="AAPL",start = dateStart, end=dateEnd, quote = c("AdjClose"),
retclass="zoo")
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
## time series starts 2005-01-03
## time series ends 2018-11-27
# Microsoft data
MSFT <- get.hist.quote(instrument="msft",start = dateStart,end=dateEnd,quote = c("AdjClose"),
retclass="zoo")
## time series starts 2005-01-03
## time series ends 2018-11-27
# IBM data
IBM <- get.hist.quote(instrument="ibm",start = dateStart,end=dateEnd,quote = c("AdjClose"),
retclass="zoo")
## time series starts 2005-01-03
## time series ends 2018-11-27
# S&P 500 index, which is one of the leading stock market indices for US equity
SP500 <- get.hist.quote(instrument="^gspc",start = dateStart, end=dateEnd,quote = c("AdjClose"),
retclass="zoo")
## time series starts 2005-01-03
## time series ends 2018-11-27
#Google data
GOOG <- get.hist.quote(instrument="GOOG",start = dateStart, end=dateEnd,quote = c("AdjClose"),
retclass="zoo")
## time series starts 2005-01-03
## time series ends 2018-11-27
data <- merge(SP500,IBM,MSFT,APPLE,GOOG) # closing price data
return.cc = diff(log(data))
ret.cc<-as.data.frame(tail(return.cc, 3000)) # latest 3000 observations
p <- seq(0.001, 0.01, 0.0001) #p for VaR
main.names<-c("SP500", "IBM", "MSFT","AAPL","GOOG")
#sample correlation
rho.cal<-function(X){
rho.hat<-cor(sign(X-mean(X)), X-mean(X))
return(rho.hat)
}
ret.0.cc<-ret.cc
for(j in 1:5){
ret.0.cc[, j]<-ret.cc[, j]-mean(ret.cc[, j])
}
rho<-apply(as.matrix(ret.0.cc), MARGIN=2, FUN=rho.cal)
#calculate degree of freedom
nu<-rep(0, 5)
for(i in 1:5){
fun <- function (x) rho[i]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
nu[i] <- uniroot(fun, c(2, 8))$root
}
number<-ncol(data)
acf.s<-rep(0, number)
acf.abs<-rep(0, number)
acf.sq<-rep(0, number)
for(j in 1:number){
acf.s[j]<-acf(ret.0.cc[, j], plot=FALSE)$acf[2]
acf.abs[j]<-acf(abs(ret.0.cc[, j]), plot=FALSE)$acf[2]
acf.sq[j]<-acf(ret.0.cc[, j]^2, plot=FALSE)$acf[2]
}
corr<-data.frame(apply(ret.cc, 2, mean), apply(ret.cc, 2, sd),
apply(ret.cc, 2, kurtosis),apply(ret.cc, 2, skewness),acf.s, acf.abs,
acf.sq, rho, nu)
rownames(corr)<-main.names
colnames(corr)<-c("mean", "sd","kurtosis","skewness","series", "abs", "sq", "sign-rho", "df")
xtable(corr, digits=4)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:18 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrrr}
## \hline
## & mean & sd & kurtosis & skewness & series & abs & sq & sign-rho & df \\
## \hline
## SP500 & 0.0002 & 0.0124 & 10.8749 & -0.3711 & -0.1002 & 0.2861 & 0.2152 & 0.6351 & 2.9877 \\
## IBM & 0.0002 & 0.0139 & 5.9892 & -0.2393 & -0.0075 & 0.2125 & 0.1483 & 0.6894 & 3.6302 \\
## MSFT & 0.0005 & 0.0172 & 9.4608 & 0.1642 & -0.0687 & 0.2374 & 0.1571 & 0.6721 & 3.3673 \\
## AAPL & 0.0009 & 0.0199 & 7.2219 & -0.4408 & 0.0077 & 0.2487 & 0.1832 & 0.6920 & 3.6775 \\
## GOOG & 0.0005 & 0.0180 & 11.6231 & 0.5192 & -0.0115 & 0.1730 & 0.0734 & 0.6622 & 3.2461 \\
## \hline
## \end{tabular}
## \end{table}
ret.0<-tail(ret.0.cc, 1000) #latest 1000 observations
##### normal garch
VaR_rateG=matrix(0, nrow=length(p), ncol=ncol(ret.0))
ES_rateG=matrix(0, nrow=length(p), ncol=ncol(ret.0))
sigma_rate=c()
omega=c()
alpha=c()
beta=c()
for (j in 1:ncol(ret.0))
{
g = garchFit(~garch(1,1),ret.0[,j],cond.dist="norm",include.mean=FALSE,
trace=FALSE)
omega[j] = g@fit$matcoef[1,1]
alpha[j] = g@fit$matcoef[2,1]
beta[j] = g@fit$matcoef[3,1]
sigma_rate[j] = omega[j] + alpha[j] * ret.0.cc[,j][T]^2 + beta[j] * g@h.t[T]
for(i in 1:length(p)){
VaR_rateG[i,j] = -sqrt(sigma_rate[j]) * qnorm(p[i]) * 1000
ES_rateG[i,j] = -sqrt(sigma_rate[j])*integrate(function(q){q*dnorm(q)},-Inf,qnorm(p[i]))$value/(p[i]) * 1000
}
}
####### t garch
VaR_ratetG1=matrix(0, nrow=length(p), ncol=ncol(ret.0))
sigma_rate=c()
omega=c()
alphaG=c()
beta=c()
ES_ratetG1=matrix(0, nrow=length(p), ncol=ncol(ret.0))
dfG=c()
for (j in 1:ncol(ret.0))
{
g = garchFit(~garch(1,1),ret.0[,j],cond.dist="std",include.mean=FALSE,
trace=FALSE)
omega[j] = g@fit$matcoef[1,1]
alphaG[j] = g@fit$matcoef[2,1]
beta[j] = g@fit$matcoef[3,1]
dfG[j]=g@fit$matcoef[4,1]
sigma_rate[j] = omega[j] + alpha[j] * ret.0[,j][T]^2 + beta[j] * g@h.t[T]
for(i in 1:length(p)){
VaR_ratetG1[i,j] =(-1)* qstd (p[i], mean = 0, sd = sqrt(sigma_rate[j]), nu = dfG[j])*1000
# using fGarch package to get df
ES_ratetG1[i,j] = sqrt(sigma_rate[j])*sqrt((dfG[j]-2)/dfG[j])*dt(qt (p[i],dfG[j]), dfG[j])/p[i]*(dfG[j] + (qt (p[i],dfG[j]))^2)/(dfG[j]-1)*1000}
}
ret.0.cc<-as.data.frame(ret.cc)
for(j in 1:5){
ret.0.cc[, j]<-ret.cc[, j]-mean(ret.cc[,j])
}
alpha.choose<-function(ret.cen, alpha=seq(0.01, 0.3, 0.01), cut.t){
rho<-rho.cal(ret.cen)
vol<-abs(ret.cen)/rho
t<-length(ret.cen)
MSE_alpha<-rep(0, length(alpha))
for(a in 1:length(alpha)){
s<-mean(vol[1:cut.t])
error<-rep(0, t)
for(i in 1:t){
error[i]<-vol[i]-s
s<-alpha[a]*vol[i]+(1-alpha[a])*s
}
MSE_alpha[a]<-mean(error[-(1:cut.t)]^2)
}
rmse<-sqrt(min(MSE_alpha))
alpha.opt<-alpha[which.min(MSE_alpha)]
return(c(rmse, alpha.opt))
}
alpha.opt<-rep(0, number)
rmse.usual<-rep(0, number)
for(num in 1:number){
result<-alpha.choose(ret.0.cc[, num], seq(0.01, 0.3, 0.01), 2000)
alpha.opt[num]<-result[2]
rmse.usual[num]<-result[1]
}
alpha.opt;rmse.usual
## [1] 0.14 0.05 0.07 0.07 0.06
## [1] 0.008803848 0.013319630 0.015985200 0.014808302 0.015813455
rmse.elastic.cal<-function(ret.cen, alpha, omega, lambda, cut.t){
rho<-rho.cal(ret.cen)
vol<-abs(ret.cen)/rho
#vol.sub<-vol-sd(ret.cen)
s<-mean(vol[1:cut.t])
sd.r<-sd(ret.cen)
if(abs(s-sd.r)>=omega*lambda)
s.elastic<- sign(s-sd.r)*(abs(s-sd.r)-omega*lambda)/(1+(1-omega)*lambda)+sd.r
else s.elastic<-sd.r
t<-length(ret.cen)
error<-rep(0, t)
for(i in 1:t){
error[i]<-vol[i]-s.elastic
s<-alpha*vol[i]+(1-alpha)*s
if(abs(s-sd.r)>=omega*lambda)
s.elastic<-sign(s-sd.r)*(abs(s-sd.r)-omega*lambda)/(1+(1-omega)*lambda)+sd.r
else s.elastic<-sd.r
}
rmse<-sqrt(mean(error[-(1:cut.t)]^2))
return(rmse)
}
lambda<-seq(0, 0.004, 0.0002)
omega<-seq(0, 1, 0.1)
rmse.elastic.opt<-rep(0, number)
omega.opt<-rep(0, number)
lambda.opt<-rep(0, number)
for(num in 1:number){
rmse.elastic<-matrix(0, nrow=length(omega), ncol=length(lambda))
for(a in 1:length(omega)){
for(l in 1:length(lambda)){
rmse.elastic[a, l]<-rmse.elastic.cal(ret.0.cc[, num], alpha.opt[num],
omega[a], lambda[l], 2000)
}
}
index.o<-which.min(apply(rmse.elastic, 1, min))
index.l<-which.min(apply(rmse.elastic, 2, min))
omega.opt[num]<-omega[index.o]
lambda.opt[num]<-lambda[index.l]
rmse.elastic.opt[num]<-rmse.elastic[index.o, index.l]
}
omega.opt; lambda.opt
## [1] 0.2 0.2 0.6 0.3 0.3
## [1] 0.0034 0.0040 0.0022 0.0034 0.0036
rmse.elastic.opt
## [1] 0.008775514 0.013299324 0.015928461 0.014782626 0.015784693
vol.forecast<-function(ret.cen, alpha, omega, lambda, cut.t){
rho<-rho.cal(ret.cen)
vol<-abs(ret.cen)/rho
sd.r<-sd(ret.cen)
t<-length(ret.cen)
res<-rep(0, t)
res.elastic<-rep(0, t)
s<-mean(vol[(1:cut.t)])
s.elastic<-sign(s-sd.r)*max((abs(s-sd.r)-
omega*lambda), 0)/(1+(1-omega)*lambda)+sd.r
fore<-rep(0, t)
fore.elastic<-rep(0, t)
for(i in 1:t){
res[i]<-ret.cen[i]/s
res.elastic[i]<-ret.cen[i]/s.elastic
fore[i]<-s
fore.elastic[i]<-s.elastic
s<-alpha*vol[i]+(1-alpha)*s
s.elastic<-sign(s-sd.r)*
max((abs(s-sd.r)-omega*lambda), 0)/(1+(1-omega)*lambda)+sd.r
}
fore.vol<-s
fore.vol.elastic<-s.elastic
return(list(fore.vol, fore.vol.elastic, res, res.elastic, fore, fore.elastic))
}
forecast.vol<-rep(0, number)
forecast.vol.elastic<-rep(0, number)
res<-ret.0.cc
res.elastic<-ret.0.cc
fore.vary<-ret.0.cc
fore.elastic.vary<-ret.0.cc
for(num in 1:number){
result<-vol.forecast(ret.0.cc[, num], alpha.opt[num], omega.opt[num],
lambda.opt[num], 2000)
forecast.vol[num]<-result[[1]]
forecast.vol.elastic[num]<-result[[2]]
res[, num]<-result[[3]]
res.elastic[, num]<-result[[4]]
fore.vary[, num]<-result[[5]]
fore.elastic.vary[, num]<-result[[6]]
}
forecast.vol; forecast.vol.elastic
## [1] 0.01462323 0.01935135 0.02464491 0.03018339 0.02447984
## [1] 0.01393917 0.01853662 0.02331948 0.02914136 0.02338628
vol.obs<-ret.0.cc
for(num in 1:number){
vol.obs[, num]<-abs(ret.0.cc[, num])/rho[num]
}
for(num in 1:number){
plot(1:1000, tail(vol.obs[, num], 1000), xlab="t", ylab="vol", main=main.names[num], lwd=3,
col="yellow", type="l")
lines(1:1000, tail(fore.vary[, num], 1000), col="blue", lwd=5)
lines(1:1000, tail(fore.elastic.vary[, num], 1000), col="purple")
abline(h=sd(ret.0.cc[, num]), lty=3, col="red", lwd=3)
legend("top",legend=c("vol", "NP", "Elastic", "sd"), col=c("yellow", "blue", "purple", "red"),
lty=c(1, 1, 1, 3), lwd=c(3, 5, 1, 3), cex=0.7)
}
for(num in 1:number){
plot(1:1000, tail(fore.vary[, num], 1000), col="green", lwd=5, type="l")
lines(1:1000, tail(fore.elastic.vary[, num], 1000), col="purple", lwd=4)
abline(h=sd(ret.0.cc[, num]), lty=3, col="red", lwd=3)
legend("top",legend=c( "NP", "Elastic", "sd"), col=c( "green", "purple", "red"),
lty=c(1, 1, 3), lwd=c( 5, 1, 3), cex=0.7)
}
vol.CI<-matrix(0, nrow=number, ncol=2)
vol.CI.elastic<-matrix(0, nrow=number, ncol=2)
for(i in 1:number){
vol.CI[i, ]<-c(forecast.vol[i]+qnorm(0.05, 0, 1)*rmse.usual[i]/sqrt(1000),
forecast.vol[i]+qnorm(0.95, 0, 1)*rmse.usual[i]/sqrt(1000))
vol.CI.elastic[i, ]<-c(forecast.vol.elastic[i]+qnorm(0.05, 0,1)*rmse.elastic[i]/sqrt(1000),
forecast.vol.elastic[i]+qnorm(0.95, 0, 1)*rmse.elastic[i]/sqrt(1000))
}
scale.p<-function(nu, p) sqrt((nu-2)/nu)*dt(qt (p,nu), nu)/p*(nu + (qt (p,nu))^2)/(nu-1)
VaR.dd<-matrix(0, nrow=length(p), ncol=number)
ES.dd<-matrix(0, nrow=length(p), ncol=number)
VaR.dd.elastic<-matrix(0, nrow=length(p), ncol=number)
ES.dd.elastic<-matrix(0, nrow=length(p), ncol=number)
for(num in 1:number){
for(i in 1:length(p)){
VaR.dd[i, num]<-forecast.vol[num]*qstd(p=p[i], nu=nu[num])*(-1000)
VaR.dd.elastic[i, num]<-forecast.vol.elastic[num]*qstd(p=0.01, nu=nu[num])*(-1000)
ES.dd[i, num]<-forecast.vol[num]*scale.p(nu[num], p[i])*1000
ES.dd.elastic[i, num]<-forecast.vol.elastic[num]*scale.p(nu[num], p[i])*1000
}
}
VaR.dd[91,]; VaR.dd.elastic[91, ]
## [1] 38.30383 51.44188 65.44769 80.22095 64.86404
## [1] 36.51201 49.27607 61.92785 77.45145 61.96643
ES.dd[91, ]
## [1] 59.18422 73.77755 96.32483 114.57789 96.79783
no.penalty<-data.frame(alpha.opt, rmse.usual, forecast.vol,(t(VaR.dd)[,91]), (t(ES.dd)[,91]) )
colnames(no.penalty)<-c("alpha.opt", "RMSE", "Volatility", "VaR", "ES")
rownames(no.penalty)<-main.names
xtable(no.penalty, digits=6)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & alpha.opt & RMSE & Volatility & VaR & ES \\
## \hline
## SP500 & 0.140000 & 0.008804 & 0.014623 & 38.303826 & 59.184225 \\
## IBM & 0.050000 & 0.013320 & 0.019351 & 51.441877 & 73.777549 \\
## MSFT & 0.070000 & 0.015985 & 0.024645 & 65.447693 & 96.324834 \\
## AAPL & 0.070000 & 0.014808 & 0.030183 & 80.220948 & 114.577890 \\
## GOOG & 0.060000 & 0.015813 & 0.024480 & 64.864037 & 96.797826 \\
## \hline
## \end{tabular}
## \end{table}
penalty<-data.frame(omega.opt, lambda.opt,rmse.elastic.opt, forecast.vol.elastic,
t(VaR.dd.elastic)[, 91],(t(ES.dd.elastic)[,91]))
colnames(penalty)<-c("omega.opt", "lambda.opt","RMSE", "Volatility", "VaR", "ES")
rownames(penalty)<-main.names
xtable(penalty, digits=6)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrr}
## \hline
## & omega.opt & lambda.opt & RMSE & Volatility & VaR & ES \\
## \hline
## SP500 & 0.200000 & 0.003400 & 0.008776 & 0.013939 & 36.512009 & 56.415641 \\
## IBM & 0.200000 & 0.004000 & 0.013299 & 0.018537 & 49.276072 & 70.671367 \\
## MSFT & 0.600000 & 0.002200 & 0.015928 & 0.023319 & 61.927854 & 91.144393 \\
## AAPL & 0.300000 & 0.003400 & 0.014783 & 0.029141 & 77.451450 & 110.622275 \\
## GOOG & 0.300000 & 0.003600 & 0.015785 & 0.023386 & 61.966433 & 92.473676 \\
## \hline
## \end{tabular}
## \end{table}
VaR.dd.CI.m<-matrix(0, nrow=number, ncol=4)
ES.dd.CI.m<-matrix(0, nrow=number, ncol=4)
for(num in 1:number){
VaR.dd.CI.m[num, 1:2]<-(-1000)*qstd(p=0.01, nu=nu[num])*vol.CI[num,]
VaR.dd.CI.m[num, 3:4]<-(-1000)*qstd(p=0.01, nu=nu[num])*vol.CI.elastic[num, ]
ES.dd.CI.m[num, 1:2]<-(1000)*scale.p(nu[num], 0.01)*vol.CI[num, ]
ES.dd.CI.m[num, 3:4]<-(1000)*scale.p(nu[num], 0.01)*vol.CI.elastic[num, ]
}
rho.res<-apply(res, 2, rho.cal)
rho.res.elastic<-apply(res.elastic, 2, rho.cal)
nu.res<-rho.res
nu.res.elastic<-rho.res.elastic
for(num in 1:number){
fun <- function (x) rho.res[num]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
nu.res[num] <- uniroot(fun, c(2, 8))$root
fun <- function (x) rho.res.elastic[num]*(x-1)*beta(x/2,1/2)-2*sqrt(x-2)
nu.res.elastic[num] <- uniroot(fun, c(2, 8))$root
}
nu.res.elastic
## Adjusted.SP500 Adjusted.IBM Adjusted.MSFT Adjusted.APPLE Adjusted.GOOG
## 4.649606 3.927487 4.045986 4.644023 3.728411
nu.res
## Adjusted.SP500 Adjusted.IBM Adjusted.MSFT Adjusted.APPLE Adjusted.GOOG
## 4.323821 3.765093 3.859418 4.472556 3.652099
VaR.dd.res<-matrix(0, nrow=length(p), ncol=number)
ES.dd.res<-matrix(0, nrow=length(p), ncol=number)
VaR.dd.elastic.res<-matrix(0, nrow=length(p), ncol=number)
ES.dd.elastic.res<-matrix(0, nrow=length(p), ncol=number)
for(num in 1:number){
for(i in 1:length(p)){
VaR.dd.res[i, num]<-forecast.vol[num]*qstd(p=p[i], nu=nu.res[num])*(-1000)
VaR.dd.elastic.res[i, num]<-forecast.vol.elastic[num]*qstd(p=p[i], nu=nu.res.elastic[num])*(-1000)
ES.dd.res[i, num]<-forecast.vol[num]*scale.p(nu.res[num], p[i])*1000
ES.dd.elastic.res[i, num]<-forecast.vol.elastic[num]*scale.p(nu.res.elastic[num], p[i])*1000
}
}
nu.res
## Adjusted.SP500 Adjusted.IBM Adjusted.MSFT Adjusted.APPLE Adjusted.GOOG
## 4.323821 3.765093 3.859418 4.472556 3.652099
VaR.dd.elastic.res[91,]
## [1] 36.55100 49.15726 61.74651 76.42124 62.13669
ES.dd.elastic.res[91,]
## [1] 49.08741 68.84120 85.76310 102.65919 88.36805
VaR.dd.CI.res<-matrix(0, nrow=number, ncol=4)
ES.dd.CI.res<-matrix(0, nrow=number, ncol=4)
for(num in 1:number){
VaR.dd.CI.res[num, 1:2]<-(-1000)*qstd(p=0.01, nu=nu.res[num])*vol.CI[num,]
VaR.dd.CI.res[num, 3:4]<-(-1000)*qstd(p=0.01, nu=nu.res.elastic[num])*vol.CI.elastic[num, ]
ES.dd.CI.res[num, 1:2]<-(1000)*scale.p(nu.res[num], 0.01)*vol.CI[num, ]
ES.dd.CI.res[num, 3:4]<-(1000)*scale.p(nu.res.elastic[num], 0.01)*vol.CI.elastic[num, ]
}
VaR.CI<-data.frame(VaR.dd.CI.m[, 1:2], VaR.dd.CI.res[, 1:2],
VaR.dd.CI.res[, 3:4], VaR.dd.CI.res[, 3:4])
xtable(VaR.CI)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & X1 & X2 & X1.1 & X2.1 & X1.2 & X2.2 & X1.3 & X2.3 \\
## \hline
## 1 & 37.10 & 39.50 & 37.35 & 39.76 & 34.39 & 38.71 & 34.39 & 38.71 \\
## 2 & 49.60 & 53.28 & 49.56 & 53.24 & 46.98 & 51.34 & 46.98 & 51.34 \\
## 3 & 63.24 & 67.66 & 63.20 & 67.61 & 59.57 & 63.92 & 59.57 & 63.92 \\
## 4 & 78.17 & 82.27 & 77.36 & 81.41 & 74.26 & 78.58 & 74.26 & 78.58 \\
## 5 & 62.68 & 67.04 & 62.88 & 67.26 & 59.95 & 64.32 & 59.95 & 64.32 \\
## \hline
## \end{tabular}
## \end{table}
ES.CI<-data.frame(ES.dd.CI.m[, 1:2], ES.dd.CI.res[, 1:2],
ES.dd.CI.res[, 3:4], ES.dd.CI.res[, 3:4])
rownames(ES.CI)<-main.names
xtable(ES.CI)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & X1 & X2 & X1.1 & X2.1 & X1.2 & X2.2 & X1.3 & X2.3 \\
## \hline
## SP500 & 57.33 & 61.04 & 51.00 & 54.30 & 46.19 & 51.98 & 46.19 & 51.98 \\
## IBM & 71.14 & 76.42 & 70.27 & 75.49 & 65.79 & 71.90 & 65.79 & 71.90 \\
## MSFT & 93.08 & 99.57 & 88.95 & 95.16 & 82.74 & 88.79 & 82.74 & 88.79 \\
## AAPL & 111.65 & 117.50 & 104.80 & 110.29 & 99.76 & 105.56 & 99.76 & 105.56 \\
## GOOG & 93.55 & 100.05 & 90.01 & 96.27 & 85.26 & 91.48 & 85.26 & 91.48 \\
## \hline
## \end{tabular}
## \end{table}
table.usual<-data.frame(alpha.opt, rmse.usual, forecast.vol, VaR.dd[length(p),], ES.dd[length(p),])
table.elastic<-data.frame(omega.opt, rmse.elastic.opt, forecast.vol.elastic, VaR.dd.elastic[length(p), ], ES.dd.elastic[length(p),])
model.risk.VaR <-matrix(0, ncol=number, nrow=length(p))
model.risk.ES <-matrix(0, ncol=number, nrow=length(p))
for(num in 1:number){
VaR<-data.frame(VaR_rateG[, num], VaR_ratetG1[, num], VaR.dd[, num], VaR.dd.res[, num])
ES<-data.frame(ES_rateG[, num], ES_ratetG1[, num], ES.dd[, num], ES.dd.res[, num])
model.risk.VaR[, num]<-apply(VaR, 1, max)/apply(VaR, 1, min)
model.risk.ES[, num]<-apply(ES, 1, max)/apply(ES, 1, min)
}
model.risk.VaR.elastic<-matrix(0, ncol=number, nrow=length(p))
model.risk.ES.elastic <-matrix(0, ncol=number, nrow=length(p))
for(num in 1:number){
VaR<-data.frame(VaR_rateG[, num], VaR_ratetG1[, num], VaR.dd.elastic[, num], VaR.dd.elastic.res[, num])
ES<-data.frame(ES_rateG[, num], ES_ratetG1[, num], ES.dd.elastic[, num], ES.dd.elastic.res[, num])
model.risk.VaR.elastic[, num]<-apply(VaR, 1, max)/apply(VaR, 1, min)
model.risk.ES.elastic[, num]<-apply(ES, 1, max)/apply(ES, 1, min)
}
MR.usual<-data.frame(model.risk.VaR[length(p),], model.risk.ES[length(p),])
colnames(MR.usual)<-c("VaR", "ES")
rownames(MR.usual)<-main.names
MR.usual
## VaR ES
## SP500 2.186598 2.929684
## IBM 1.831004 2.292130
## MSFT 2.078809 2.670553
## AAPL 2.378313 2.964999
## GOOG 1.931144 2.507516
MR.elastic<-data.frame(model.risk.VaR.elastic[length(p),], model.risk.ES.elastic[ length(p),])
colnames(MR.elastic)<-c("VaR", "ES")
rownames(MR.elastic)<-main.names
MR.elastic
## VaR ES
## SP500 2.072867 2.792636
## IBM 1.753915 2.195627
## MSFT 1.967009 2.526928
## AAPL 2.296205 2.862637
## GOOG 1.844097 2.395500
MR<-data.frame(MR.usual, MR.elastic)
xtable(MR)
## % latex table generated in R 4.0.0 by xtable 1.8-4 package
## % Mon Jul 27 15:24:24 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & VaR & ES & VaR.1 & ES.1 \\
## \hline
## SP500 & 2.19 & 2.93 & 2.07 & 2.79 \\
## IBM & 1.83 & 2.29 & 1.75 & 2.20 \\
## MSFT & 2.08 & 2.67 & 1.97 & 2.53 \\
## AAPL & 2.38 & 2.96 & 2.30 & 2.86 \\
## GOOG & 1.93 & 2.51 & 1.84 & 2.40 \\
## \hline
## \end{tabular}
## \end{table}
temp1 <- data.frame (p, as.data.frame(model.risk.VaR))
colnames (temp1) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp1 <- gather (temp1, "Stock", "Model_risk", 2:6)
temp2 <- data.frame (p, as.data.frame(model.risk.ES))
colnames (temp2) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp2 <- gather (temp2, "Stock", "Model_risk", 2:6)
plot1 <- ggplot(temp1) +
geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
ylab ("Model risk of VaR: no penalty")+ylim(c(1.5, 4))
plot2 <- ggplot(temp2) +
geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
ylab ("Model risk of ES: no penalty")+ylim(c(1.5, 5))
grid.arrange(plot1, plot2, ncol=2)
temp1 <- data.frame (p,as.data.frame(model.risk.VaR.elastic))
colnames (temp1) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp1 <- gather (temp1, "Stock", "Model_risk", 2:6)
temp2 <- data.frame (p, as.data.frame(model.risk.ES.elastic))
colnames (temp2) <- c ("p","S&P 500","IBM","MSFT","AAPL", "GOOG")
temp2 <- gather (temp2, "Stock", "Model_risk", 2:6)
plot1 <- ggplot(temp1) +
geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
ylab ("Model risk of VaR: elastic net")+ylim(c(1.5, 4))
plot2 <- ggplot(temp2) +
geom_point (aes (x = p, y = Model_risk, color = Stock, shape = Stock)) +
ylab ("Model risk of ES: elastic net")+ylim(c(1.5, 5))
grid.arrange(plot1, plot2, ncol=2)
model.risk.V<-list()
par(mfrow=c(1, 5))
for(j in 1:5){
model.risk.V[[j]]<-data.frame(model.risk.VaR.elastic[,j], model.risk.VaR[, j])
colnames(model.risk.V[[j]])<-c("REG", "NR")
boxplot(model.risk.V[[j]], col=c(2, 3), main=main.names[j])
}
model.risk.E<-list()
par(mfrow=c(1, 5))
for(j in 1:5){
model.risk.E[[j]]<-data.frame(model.risk.ES.elastic[,j ], model.risk.ES[, j])
colnames(model.risk.E[[j]])<-c("REG", "NG")
boxplot(model.risk.E[[j]], col=c("yellow", 5), main=main.names[j])
}
model.risk.0.01.nonPen<-data.frame(model.risk.VaR[length(p),], model.risk.ES[ length(p),])
model.risk.0.01<-data.frame(model.risk.VaR.elastic[length(p),], model.risk.ES.elastic[ length(p),], model.risk.0.01.nonPen)
rownames(model.risk.0.01)<-main.names
colnames(model.risk.0.01)<-c("VaR(REG)", "ES(REG)", "VaR", "ES")
#library(vioplot)
par(mfrow=c(1, 5))
for(i in 1:number){
vioplot(model.risk.V[[i]], col=c(2, 3), names=c("REG", "NR"))
title(paste(main.names[i], "(VaR)"), ylab="model risk")
}
par(mfrow=c(1, 5))
for(i in 1:number){
vioplot(model.risk.E[[i]], col=c("yellow", 5), names=c("REG", "NR"))
title(paste(main.names[i], "(ES)"), ylab="model risk")
}
EVT for VaR and ES forecasting: Reg:
res<-tail(res, 1000)
res.elastic<-tail(res.elastic, 1000)
m<-seq(150, 300, 10)
m.opt<-rep(0, number)
for(j in 1:number){
res.use<-res[, num]
mse<-rep(0, length(m))
for(i in 1:length(m)){
res.order<-sort(res.use)[1:m[i]]
y<-log((-1)*res.order)
x<-log(seq(1, m[i], 1)/1000)
fit<-lm(y~x)
mse[i]<-mean(fit$residuals^2)
}
m.opt[j]<-m[which.min(mse)]
}
m.opt
## [1] 150 150 150 150 150
m.opt.elastic<-rep(0, number)
for(j in 1:number){
res.use<-res.elastic[, num]
mse<-rep(0, length(m))
for(i in 1:length(m)){
res.order<-sort(res.use)[1:m[i]]
y<-log((-1)*res.order)
x<-log(seq(1, m[i], 1)/1000)
fit<-lm(y~x)
mse[i]<-mean(fit$residuals^2)
}
m.opt.elastic[j]<-m[which.min(mse)]
}
m.opt.elastic
## [1] 150 150 150 150 150
a.reg.inv<-rep(0, number)
a.reg.inv.elastic<-rep(0, number)
for(j in 1:number){
res.use<-sort(res[, j])[1:m.opt[j]]
y<-log((-1)*res.use)
x<-log(seq(1, m.opt[j], 1)/1000)
fit<-lm(y~x)
a.reg.inv[j]<-(-1)*fit$coefficients[2]
}
for(j in 1:number){
res.use<-sort(res.elastic[, j])[1:m.opt.elastic[j]]
y<-log((-1)*res.use)
x<-log(seq(1, m.opt.elastic[j], 1)/1000)
fit<-lm(y~x)
a.reg.inv.elastic[j]<-(-1)*fit$coefficients[2]
}
a.reg.inv
## [1] 0.4856916 0.4728223 0.4574394 0.3882623 0.4347708
a.reg.inv.elastic
## [1] 0.4673441 0.4725276 0.4498897 0.3823091 0.4295732
p<-0.1
VaR.reg<-rep(0, number)
VaR.reg.elastic<-rep(0, number)
for(j in 1:number){
VaR.reg[j]<-(-1000)*sort(res[, j])[100]*forecast.vol[j]*(0.1/0.01)^a.reg.inv[j]
VaR.reg.elastic[j]<-(-1000)*sort(res.elastic[, j])[100]*
forecast.vol.elastic[j]*(0.1/0.01)^a.reg.inv.elastic[j]
}
VaR.reg; VaR.reg.elastic
## [1] 47.24521 65.77080 75.12274 83.91127 69.36910
## [1] 42.80098 59.26378 65.42395 76.52281 62.27243
ES.reg<-VaR.reg*(1/(1-a.reg.inv))
ES.reg.elastic<-VaR.reg.elastic*(1/(1-a.reg.inv.elastic))
Hill:
nc<-seq(5, 300, 5)
hill.est<-function(y, nc){
y.use<-sort(y)[1:nc]
a.hill<-nc/sum(log(y.use/y.use[nc]))
return(a.hill)
}
for(j in 1:number){
a.hill<-rep(0, length(nc))
for(i in 1:length(nc)){
a.hill[i]<-hill.est(res[, j], nc[i])
}
plot(nc, a.hill, type="b", main=main.names[j])
}
nc.opt<-100
for(j in 1:number){
a.hill<-rep(0, length(nc))
for(i in 1:length(nc)){
a.hill[i]<-hill.est(res.elastic[, j], nc[i])
}
plot(nc, a.hill, type="b", main=main.names[j])
}
nc.opt.elastic<-100
a.hill.inv<-rep(0, number)
a.hill.inv.elastic<-rep(0, number)
for(j in 1:number){
a.hill.inv[j]<-1/hill.est(res[, j], nc.opt)
a.hill.inv.elastic[j]<-1/hill.est(res[, j], nc.opt.elastic)
}
VaR.hill<-rep(0, number)
VaR.hill.elastic<-rep(0, number)
for(j in 1:number){
VaR.hill[j]<-(-1000)*sort(res[, j])[100]*forecast.vol[j]*(0.1/0.01)^a.hill.inv[j]
VaR.hill.elastic[j]<-(-1000)*sort(res.elastic[, j])[100]*
forecast.vol.elastic[j]*(0.1/0.01)^a.hill.inv.elastic[j]
}
ES.hill<-VaR.hill*(1/(1-a.hill.inv))
ES.hill.elastic<-VaR.hill.elastic*(1/(1-a.hill.inv.elastic))
ES.hill; ES.hill.elastic
## [1] 94.32527 101.66041 126.63569 166.99158 136.43301
## [1] 89.13979 91.66483 112.22022 154.38967 123.95010
EVT.NP<-data.frame(VaR.hill, VaR.reg, ES.hill, ES.reg)
EVT.REG<-data.frame(VaR.hill.elastic, VaR.reg.elastic, ES.hill.elastic, ES.reg.elastic)
EVT.NP; EVT.REG
## VaR.hill VaR.reg ES.hill ES.reg
## 1 47.92611 47.24521 94.32527 91.86162
## 2 58.65128 65.77080 101.66041 124.76021
## 3 71.45808 75.12274 126.63569 138.45962
## 4 93.95547 83.91127 166.99158 137.16870
## 5 73.60442 69.36910 136.43301 122.72739
## VaR.hill.elastic VaR.reg.elastic ES.hill.elastic ES.reg.elastic
## 1 45.29139 42.80098 89.13979 80.3539
## 2 52.88450 59.26378 91.66483 112.3543
## 3 63.32371 65.42395 112.22022 118.9288
## 4 86.86519 76.52281 154.38967 123.8853
## 5 66.87000 62.27243 123.95010 109.1681
EVT<-data.frame(EVT.NP, EVT.REG)
rownames(EVT)<-main.names
EVT
## VaR.hill VaR.reg ES.hill ES.reg VaR.hill.elastic VaR.reg.elastic
## SP500 47.92611 47.24521 94.32527 91.86162 45.29139 42.80098
## IBM 58.65128 65.77080 101.66041 124.76021 52.88450 59.26378
## MSFT 71.45808 75.12274 126.63569 138.45962 63.32371 65.42395
## AAPL 93.95547 83.91127 166.99158 137.16870 86.86519 76.52281
## GOOG 73.60442 69.36910 136.43301 122.72739 66.87000 62.27243
## ES.hill.elastic ES.reg.elastic
## SP500 89.13979 80.3539
## IBM 91.66483 112.3543
## MSFT 112.22022 118.9288
## AAPL 154.38967 123.8853
## GOOG 123.95010 109.1681